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Abstract 

The objective of this study has been to create a Stirling engine model for studying the effects of regenerator 
matrix temperature oscillations on Stirling engine performance. 

A one-dimensional model with axial discretisation of engine components has been formulated using the control 
volume method. The model contains a system of ordinary differential equations (ODEs) derived from mass and 
energy balances for gas filled control volumes and energy balances for regenerator matrix control masses. 
Interpolation methods with filtering properties are used for state variables at control volume interfaces to reduce 
numerical diffusion and/or non-physical oscillations. Loss mechanisms are included directly in the governing 
equations as terms in the mass and energy balances. 

Steady state periodic solutions that satisfy cyclic boundary conditions and integral conditions are calculated 
using a custom built shooting method. 

It has been found possible to accurately solve the stiff ODE system that describes the coupled thermodynamics 
of the gas and the regenerator matrix and to reliably find periodic steady state solutions to the model. Preliminary 
results indicate that the regenerator matrix temperature oscillations do have significant impact on the regenerator 
loss, the cycle power output, and the cycle efficiency and thus deserve further study. 

© 2005 Elsevier Ltd. All rights reserved. 


1. Introduction 

Computer simulation has been used for design and optimisation of Stirling engines for several 
decades. Early modelling efforts resulted in very simple models that could be solved on the computers of 
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Nomenclature 

A wetted surface area [m 2 ]. 
c p spec, heat for const, pressure [J/(kg K)] 

c y spec, heat for const, volume [J/(kg K)] 

E energy [J] 

h conv. heat transfer coeff. [W/(m 2 K)] 

m mass [kg] 

lit mass flow [kg/s] 

p pressure [Pa] 

t time [s] 

T temperature [K] 

V volume [m 3 ] 

cv in control volume 

lbnd at left boundary of control volume 

rbnd at right boundary of control volume 


the day. As computers grow ever more powerful the many assumptions in early models can be gradually 
relaxed in order to achieve a more accurate prediction of engine behaviour. With more accurate models 
new insight into phenomena occurring in engines can be gained and engine designs can be optimised to 
achieve higher thermal efficiencies. 

One classical assumption in Stirling engine models is that the temperature of all metal in the engine 
remains constant during the cycle. In one of the key components in a Stirling engine, the regenerator, the 
validity of this assumption can be questionable. Analytical studies, such as the studies by Jones [1] 
suggest that metal temperature oscillations in regenerators can have an impact on engine performance. 

Experimental studies of regenerators have mostly been conducted in order to generate correlations for 
heat transfer and flow friction, such as those compared by Thomas and Pittman [2], that can be used in 
simulation programs. Several experimental studies, for instance the work of Gedeon and Wood [3] and 
of Isshiki et al. [4], have been performed by placing small samples of regenerator matrix material into 
specialised regenerator test equipment instead of performing the measurements on complete 
regenerators inside actual engines. These studies have yielded information about the flow friction and 
heat transfer in regenerator matrices but they do not show the influence of regenerator performance on 
machine performance. In the heat transfer study by Siegel [5], measurements have been performed on a 
full regenerator inside an actual Stirling engine. But it has only been possible to measure bulk 
temperatures and pressure losses that result from the combined effects of all phenomena, such as heat 
transfer, compression/expansion, matrix temperature oscillations, and flow friction, occurring in the 
regenerator in the engine; The influence of the individual phenomena has not been revealed. Some 
authors report applying correction terms to simulation results calculated with constant regenerator 
matrix temperatures in order to account for the effects of regenerator matrix temperature oscillations, see 
for instance Jones or Kiihl and Schultz [1,6]. 

This paper presents an overview of the method and the initial findings of an effort to explore the 
effects of dynamic temperature oscillations in regenerator matrices using a numerical simulation model. 
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2. The Stirling engine 

Similarly to an internal combustion engine the Stirling engine is based on a gas cycle, where work is 
expended to compress a cold working gas and where work is extracted by expanding the gas after it has 
been heated in order to increase the pressure. The Stirling engine, however, works on a closed gas cycle, 
where the working gas does not take part in combustion and is not exchanged in every cycle. Instead the 
heating and cooling of the working gas is achieved by sending the working gas back and forth through a 
serial connection of three heat exchangers, viz. a heater, a cooler, and a regenerator, configured as 
sketched in Fig. 1. 

The efficiency of a Stirling engine can be improved by the regenerator because it can recycle 
some of the heat that is removed from the gas during transfer to the cold cylinder to preheat the 
gas when it is transferred back to the hot cylinder. The regenerator is usually made as a porous 
matrix with a large heat transfer area and a large specific heat capacity from, for instance, metal 
felt. The regenerator acts as a thermal heat storage that absorbs heat when hot gas flows through it 
and then releases it again when the flow direction is reversed and the gas flow becomes cooler than 
the matrix. At steady state periodic operating conditions there is a steep temperature gradient 
through the regenerator matrix dictated to some degree by the temperatures of the heater and 
cooler. In order to achieve a high thermal efficiency, it is very important to have a well performing 
regenerator so that the energy flux loss through the regenerator from the heater to the cooler is 
minimised. At the same time the pressure drop across the regenerator must be kept as low as 
possible in order to avoid wasting too much work on pushing the gas back and forth. Regenerator 
optimisation is thus at the heart of Stirling engine design. 

Many practical Stirling engine designs look very different from Fig. 1 but the net effects of 
compression, expansion, and pushing gas back and forth through heat exchangers and a regenerator are 
the same. 

The engine design parameters and operating conditions that are needed as input parameters for 
the Stirling engine model used in this study are based on an existing 9 kW engine design by 
Carlsen and Bovin [7]. A drawing of the region around the cylinder of this engine can be seen in 
Fig. 2. This engine uses one working piston and a displacer piston instead of the two working 
pistons in the engine of Fig. 1. 


Cold Cylinder Cooler Regenerator Heater Hot Cylinder 



Fig. 1. Stirling engine with two heat exchangers and a regenerator. 
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Heater- 

Hot cylinder volume 

Regenerator - 

Displacer piston — 

Cooler _ 

Cold cylinder volume 
Working piston - 



Fig. 2. Cylinder region of 9 kW Stirling engine. 


3. The Stirling engine model 

The Stirling engine model has been formulated as a system of ODEs using the control volume method 
and conservation of mass and energy. 

3.1. Discretisation in model 

For the sake of simplicity the geometry of the model is not referred to Fig. 2 but to the equivalent 
geometry in Fig. 1. 

The cylinder volumes and the manifold volumes between the components are represented by single 
control volumes. The sizes of the control volumes representing the cylinder volumes vary in a cyclic 
fashion. 

In the model, the regenerator and the heater are split into three sub components each. In the case of the 
heater this division has been chosen because the end sections of the heater tubes do not receive the same 
external heat input as the central part of the tubes. Treating the inactive parts of the heater as separate sub 
components simplifies the model code. The regenerator has been split into three components to facilitate 
non-uniform spatial discretisation where finer grids are used in the ends of the regenerator than in the 
central part. This helps to better resolve localised phenomena in the ends of the regenerator with minimal 
computational efforts. 

The cooler and the sub components in the regenerator and heater are further subdivided into 
parameterised numbers of control volumes. 

The control volumes in the discretisation described above are linked together into a string of control 
volumes that represents the space in the engine occupied by the working gas. 
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1 2 3 4 5 6 7 8 9 10 11 12 13 


Components: 

1: Cold cylinder volume 3: Cooler 5,6,7: Regenerator 

13: Hot cylinder volume 9, 10, 11 :Heater 2, 4, 8, 12: Manifold volumes 

Fig. 3. Discretisation used in Stirling engine model. 

A string of control masses that represents the regenerator matrix has also been defined. This string of 
control masses only spans the length of the regenerator and uses the same axial discretisation as the 
corresponding control volumes for the working gas. 

The resulting discretisation of the engine is illustrated in Fig. 3. 


3.2. Phenomena included in model 

The working gas is helium and is modelled as an ideal gas with non-constant thermophysical 
properties. 

The model has been formulated using control volumes, control masses, and conservation of mass and 
energy. Conservation of momentum has been reduced by assuming that the axial inertia of the gas is 
negligible. More accurately, the velocities of gas flows between control volumes are found using steady 
state correlations for pressure losses in tube flow, through regenerator matrices, and through flow 
constraints. 

Heat transfer is included between the gas and the wetted metal walls of all components. The 
instantaneous rates of heat exchange between gas and walls is calculated from the wall temperatures and 
the average temperatures in the gas filled control volumes using empirical correlations for heat exchange 
in tube flow, flow through porous matrices and flow in cylinder volumes. The temperatures of the walls 
in the cylinders, manifold volumes, the heater, and the cooler are assumed fixed at their cycle mean 
values. 

In the regenerator heat transfer is included between the gas and the regenerator matrix. The 
temperatures of the regenerator matrix control masses can be treated in two different ways. Either the 
temperatures can be assumed fixed at their cycle mean values or they can be modelled using ODEs 
derived from an energy balance for a lumped control mass, i.e. a control mass with a uniform 
temperature. The validity of using a lumped formulation for the matrix thread material has been verified 
using a separate model that resolves the radial temperature variations inside a thread subjected to the 
conditions in the regenerator. 

Heat conduction inside the metal walls of the engine is also included. This heat conduction affects the 
temperatures of the metal walls in contact with the gas and thus also the gas. The main heat conduction 
paths in the model are illustrated by the arrows in Fig. 4. 
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Fig. 4. Heat conduction paths in Stirling engine model. 


3.3. Model equation system 


The control volumes that represent the working gas are described by two ODEs each, viz. a mass 
balance and an energy balance. The mass balance ODEs take the form of Eq. (1). 


dm cv 

At 


^lbnd ^rbnd 



Flows across boundaries are represented as positive in the direction from left to right (from cold to hot) 
in Fig. 3. The influence of flow friction on the mass distribution in the engine is directly included in the 
mass balance equations because the mass flow rates are calculated from correlations for pressure losses. 

In the energy balances for the control volumes, it has been assumed that kinetic and gravitational 
potential energies in the control volumes are of negligible magnitudes. The total energy in the control 
volumes has thus been assumed equal to the internal energy. Because of this the energy balances for the 
control volumes take the form of Eq. (2). 


dffcv 

At 


A(m cy c y T cy ) _ \ i a rr t \ dV c\ 

^lbnd^p^lbnd ^rbnd^p^rbnd ' Ow ^cv ) P cv 


At 


At 



The first two terms in Eq. (2) represent enthalpy flows between neighbouring control volumes. Notice 
that interpolated boundary temperatures are used in these terms to avoid numerical diffusion. The 
interpolation is carried out using interpolation methods with filtering properties similar to the mixed 
interpolation and extrapolation method presented in by Kiihl and Schultz [8]. The third term represents 
convective heat exchange between gas and metal. The last term represents work from volume changes 
and is only relevant for the cylinder volumes. Energy loss terms such as metal heat conduction typically 
affect the wall temperatures in the engine and are thus included in the ODEs through the convective heat 
exchange term. 

The control masses that represent the regenerator matrix are described either by ODEs of the form 
Eq. (2) with only a heat exchange term or by a constant temperature. Heat conduction through the metal 
in the matrix is not taken into account. 
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In order to be able to find the temperatures of metal walls throughout the engine it is also necessary to 
integrate the heat transfer to the wall sections during each cycle. ODEs that have just the heat exchange 
term of Eq. (2) are used for this purpose. 

The governing equations in the model make up a coupled system of first order ODEs. 

3.4. Mathematical problem definition 

For the purposes of this investigation steady state periodic solutions to the ODE system of the model 
are needed. In a steady state periodic solution the values of all variables representing masses and 
energies must be the same at the end and the beginning of the cycle. This requirement is a cyclic 
boundary condition. The initial masses and energies in the control volumes and control masses must be 
found so that cyclic boundary conditions are satisfied. 

Cyclic heat transfer is assumed to take place for the metal walls in the cylinders, manifold volumes, 
and in the end sections of the heater. When constant regenerator matrix temperatures are assumed then 
conditions for cyclic heat transfer are also used for finding the matrix temperatures. 

Integrations of heat transfer to or from the wall sections and portions of the regenerator matrix, where 
cyclic heat transfer is assumed must result in zero net heat transfer if the prescribed constant 
temperatures are to be true steady state cycle mean values. The fixed wall and matrix temperatures must 
be found so that this is achieved. 

An infinite number of solutions satisfy the above conditions for steady state solutions for the model 
investigated here. In order to get solutions that are directly comparable, it is necessary to add an 
additional condition that fixes the total mass of the working gas in the engine. For this purpose, an 
integral condition that measures the deviation from a user specified mean pressure in the hot cylinder has 
been chosen. The amount of mass and energy in the string of control volumes representing the working 
gas must be scaled to eliminate this deviation. 

The mathematical problem that must be solved thus contains a cyclic boundary value problem for a 
system of ODEs, a set of integral conditions for cyclic heat transfer, and an integral condition that must 
be used for scaling the masses and energies in the gas filled control volumes. The numerical method must 
also allow additional integrations that can be used for integrating the heat absorbed in the heater and/or 
cooler and the work done on the pistons, so that performance characteristics, such as power output and 
thermal efficiency, can be calculated. 


4. Implementation of model and numerical method 

The mathematical problem of the model is solved numerically using the MusSim (Multi purpose 
software for Simulation) software by Andersen [9,10]. The software and the model is implemented in 
standard Fortran95 and has been tested to be compatible with Compaq, Intel, Lahey/Fujitsu, HP, IBM, 
and Sun compilers and thus runs on a large number of platforms. 

A good Stirling engine model should approach periodic steady state asymptotically if integrated 
forward in time from a given set of initial values, but the convergence may be slow and asymptotical. To 
drastically accelerate this convergence the MusSim software uses a purpose built shooting method. This 
shooting method treats the ODEs in the model as a boundary value problem that is solved for initial 
values and parameters corresponding to a periodic steady state solution to the initial value problem (IVP) 
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defined by the ODEs. The cyclic boundary conditions and the integral conditions are treated as residual 
equations that must be solved in much the same way as a Newton-Rhapson method solves non-linear 
residual equations. The main difference is that the shooting method requires some of the variables to be 
scaled based on a residual. 

The system of ODEs turns out to be very stiff. This is not surprising as the variables in the system have 
very different inertial properties. For instance, the mass of gas and the energy in a control volume where 
gas is blowing through at high velocity can change on a much smaller time scale than the thermal energy 
stored in a control mass in the regenerator matrix. In order to handle this stiffness a suitable IVP method 
must be used in the shooting method. During this investigation the semi-implicit GERK method by 
Thomsen [11] has been successfully applied. 

To be able to achieve accurate results it has proven necessary to solve ill conditioned non-linear 
equation systems with very strict tolerances in the internal iterations of the GERK method. This is 
handled using the built in non-linear equation solver of the MusSim software. It can be time consuming 
to robustly seek out difficult solutions of this kind but the MusSim software has features to speed up this 
process and this makes the model manageable for many purposes on a standard PC. 


5. Model verification 

A great deal of effort has been put into verifying the correctness of the model. 

One of the most obvious and also most important tests is to make sure that the model does not leak 
mass and/or energy due to misconnected control volumes or control masses. Testing for conservation of 
mass is trivial and has been done by monitoring the sum of the masses of gas in the control volumes. 
Testing for conservation of energy in the same manner is slightly more involved as it requires the system 
to be kept adiabatic from the surroundings and it requires the work terms from the movement of the 
pistons to be zeroed out. For steady state solutions, though, conservation of energy can be verified by an 
overall energy balance applied at the outer boundaries of the model. These tests are integrated into the 
model and are performed for every calculated periodic steady state solution. Conservation of mass and 
energy has been found to satisfy relative tolerances that are stricter than the tolerance enforced by the 
GERK method for the individual masses and energies. 

The correctness of the implementations of the correlations needed in the model have been tested by 
exchanging them with different correlations and checking the differences caused by the exchanges. The 
variations have been found to be within the expected ranges. 

The calculated performance of the engine according to the assumed best model has been compared to 
values measured on the actual 9 kW engine. Depending on the combination of correlations used to describe 
pressure losses and heat transfer the model overestimates the work output from the gas cycle by 8-21% and 
the efficiency of the gas cycle by 6-18%. These results were expected as the model does not account for 
certain known losses associated with the displacer piston that divides the cylinder volume of the engine. 


6. Preliminary results and discussion 

After initial verification of the model had been carried out it was studied how the spatial discretisation 
in the regenerator influences the calculated temperature profile for the regenerator matrix. During this 
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investigation the cooler was divided into 10 control volumes and the total length of the heater was 
divided into 20 control volumes. 

For fixed matrix temperatures the total length of the regenerator was divided into 10, 20, and 50 
control volumes of uniform size in subsequent simulations. The resulting matrix temperature profiles are 
shown in Fig. 5. The slightly curved shape of the calculated temperature profiles is largely independent 
of the fineness of the discretisation. 

This experiment was repeated with dynamic regenerator temperatures. In order to present the information 
from these simulations it has been chosen to plot the maximum and minimum energy based temperatures in 
the regenerator matrix control masses during a cycle. These plots are shown in Fig. 6 and it is apparent that 
the shapes of the profiles depend on the fineness of the discretisation. It also appears that the curvatures of the 
temperature profiles are largest in the ends of the regenerator. This large curvature in the ends of the matrix is 
caused by the matrix exchanging considerable amounts of energy with incoming gas flows whose 
temperatures differ significantly from the temperatures of the ends of the matrix. 

Based on these observations a new non-uniform discretisation with a larger density of control 
volumes in the ends of the regenerator was devised. In this discretisation the ends of the regenerator 
(components 5 and 7 in Fig. 3) were both given lengths corresponding to 8% of the total length of the 
regenerator while component 6 covered the remaining central part. Component 5 was divided into 10 
control volumes, component 6 into five control volumes, and the fineness of the discretisation in 
component 7, located where the largest gradients were observed in Fig. 6, was varied between 10 and 20 
control volumes. Temperature profiles showing the maximum and minimum temperatures in the matrix 
control masses during the cycle when using the new discretisation are shown in Fig. 7. The temperature 
profiles in Fig. 7 appear very similar and because of this the discretisation has not been further refined. 

Fig. 8 shows comparisons of the temperature profiles for fixed and dynamic matrix temperatures 
calculated using fine discretisations. The figure shows that the temperature profile in the regenerator 
matrix looks significantly different when the oscillations of the matrix temperatures are taken into 
account. The main differences are the temperature oscillations in the ends of matrix and a slightly flatter 
temperature gradient in the central part of the regenerator. 

In order to verify that the solutions with dynamic matrix temperatures found by the shooting method 
of the MusSim software are true steady state solutions a test procedure was devised. 

In this test procedure the initial values found by the shooting method are used as initial values for 
sequential simulations. These simulations run the engine for an additional 100 consecutive cycles from 



Fig. 5. Axial temperature profiles in the regenerator matrix calculated with fixed matrix temperatures. 
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Fig. 6. Max. and min. matrix temperature profiles calculated with dynamic matrix temperatures and uniform discretisation. 

the initial values found by the shooting method. After completing the 100 cycles the final values are 
recorded. The initial values found by the shooting method are then subtracted from these final values and 
the resulting differences are divided by the initial values. The results of this procedure we label as the 
relative drifts in the periodic steady state solution. 



Fig. 7. Max. and min. temperature profiles in the regenerator matrix calculated with dynamic matrix temperatures and non- 
uniform discretisation. 
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Fig. 8. Comparisons of temperature profiles for fixed and dynamic matrix temperatures. 

This test procedure has been completed for the above mentioned five discretisations used in the 
simulations with dynamic regenerator temperatures. The calculated relative drifts for all variables are 
shown in Fig. 9 for the case, where shooting has been carried out to a tolerance of 10 _ 5 using a relative 
tolerance of 10 in the GERK method. The figure shows that the relative drifts of all variables in all the 
tested solutions are below the tolerance enforced by the shooting method. The figure also shows that the 
signs of the largest drifts seem random in nature within components. They can thus be interpreted as 
noise from the integration of the ODEs. The variables most affected by this noise are the masses of gas in 
the regenerator control volumes and, to a lesser extent, the masses and energies in and close to the 
cylinder volumes. 

The visual inspection of the temperature profiles presented above illustrates that the temperature 
profiles appear to be shaped differently when oscillations in matrix temperatures are taken into account. 
The visual inspection does not, however, show that the temperature oscillations have an impact on the 
performance of Stirling engines. The effect on performance by the temperature oscillations is shown by 
Table 1 that contains characteristic numbers for the performance of the gas cycle in the engine for all the 
solutions presented above. The values in the table are the work output, the heat intake from the heater, 
the thermal efficiency, and the regenerator loss defined as the average net energy flux through the 
regenerator from the hot end towards the cold end of the regenerator. 



Fig. 9. Relative drifts in periodic steady state solutions. 
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Table 1 


Calculated performance of gas cycle in engine 



Work (kW) 

Heat input (kW) 

Efficiency (%) 

Reg. loss (kW) 

Fixed 10 

15.213 

31.938 

47.633 

1.519 

Fixed 20 

15.192 

31.886 

47.647 

1.519 

Fixed 50 

15.187 

31.873 

47.650 

1.518 

Dyn 10 

14.617 

30.738 

47.554 

0.543 

Dyn 20 

14.542 

30.622 

47.489 

0.528 

Dyn 50 

14.392 

30.478 

47.222 

0.496 

Dyn 10-5-10 

14.306 

30.474 

46.946 

0.470 

Dyn 10-5-20 

14.285 

30.456 

46.904 

0.467 


The numbers in Table 1 show several distinct differences in performance caused by taking dynamic 
temperature oscillations in the regenerator matrix into account. 

The net energy flux loss through the regenerator is reduced by roughly two thirds. In spite of the 
reduction of this loss term the thermal efficiency is reduced due to a significant reduction in the power 
output of the cycle. These differences will be further investigated. 

When comparing the above performance numbers for the gas cycle of the engine to the performance 
of this and other engines it is important to remember that the above performance numbers are for the gas 
cycle of the engine and not for the shaft power from the engine or for the electrical power output from an 
engine-generator assembly. 


7. Conclusion 

It has been found that true steady state periodic solutions can be reliably calculated for a cyclic 
boundary value problem that describes a Stirling engine and includes the coupled thermodynamics of a 
gas and a regenerator matrix. 

It has also been found that the calculated temperature profile in a regenerator matrix can look 
significantly different when matrix temperature oscillations are taken into account. The main differences 
have been found to be significant matrix temperature oscillations in the ends of the regenerator and a less 
steep matrix temperature gradient in the central part of the regenerator. 

Finally, it has been found that the oscillations in the regenerator matrix temperatures influence the 
calculated performance of the gas cycle in the engine. The observed effects include a reduction of the 
regenerator loss, the power output, and the thermal efficiency. 

Future work will include improvements to the model and the MusSim software as well as 
comparisons of obtained simulation results with measurements on a test bench. 
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